%% This function estimate the Chen and Szroeter (2012) test
%% Inputs:
%mu=estimeted constraints, 
%V= estimated asymptotic variance of sqrt(n)*(mu-mu_0)
%n=sample size
%P: set P=1 for step-at-unity P=2 for normal and P=3 for logistic, 
%K=tuning parameter, possible choices: sqrt(n/log(n))and sqrt(n/(2*log(log(n)))). In the paper we use the latter. 
%cl= confidence level

%% Outputs
% Rc=1 if the test rejects the null
% Pc=P-value

function [Rc,Pc]=chsztest(mu,V,n,P,K,cl)



% Create a vector which contains the invers of the sqrt of the main diagonal of the asymptotic variance
theta=1./sqrt(diag(V));
mu=mu(theta<Inf & theta>-Inf & theta ~=0);
V=V(theta<Inf & theta>-Inf & theta ~=0,theta<Inf & theta>-Inf & theta ~=0);
theta=theta(theta<Inf & theta>-Inf & theta ~=0);

switch P
    case 1 % Step-at-unity
        psi=K*theta.*mu<=1;
        % adjustment term
        Lambda =-normpdf(sqrt(n)/K)*(ones(length(mu),1));
        
    case 2 % Normal
        psi=1-normcdf(K*theta.*mu);
        % adjustment term
        Lambda =(-normpdf(K*theta.*mu)*K)/sqrt(n);
    case 3 % Logistic
        psi=(1+exp(K*theta.*mu)).^(-1);
        % adjustment term
        Lambda =((-exp(K*theta.*mu)./(1+exp(K*theta.*mu)).^2))*(K/sqrt(n));
end

% Create a diagonal matrix with the standard deviations of the elements of mu
delta=diag(theta);
% Compute the numerator of the test statistics
Q1=sqrt(n)*psi'*delta*mu-ones(1,length(mu))*Lambda;
% Compute the denominator of the test statistics
Q2=sqrt(psi'*delta*V*delta*psi);
% Compute the P-value
if Q2==0 
    Pc=1;
else
    Pc=normcdf(Q1/Q2);
end
% Reject if Pc<cl
Rc=Pc<cl;